////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand, 31 July 2013 -----------------                        //
//                                                                                                        //
//  The routine "GaussSeidel" solves the 2D Laplace equation using Gauss-Seidel relaxation method.                   //
//  Dirichlet boundary condition is used. The region is a rectangle of side lengths Lx and Ly.            //
//                                                                                                        //
//                                                                                                        //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

const Nxsize=100,Nysize=200;

int i,j,n,Nx=100,Ny=200,T=15000;

double Lx=1,Ly=2,a1=0.,a2=0.,a3=1.,a4=0.,delx,dely;

double Phi[Nxsize+1][Nysize+1],Phinew[Nxsize+1][Nysize+1]; // 2D arrays "Phi" and "Phinew" store the current and updated values of potential at grid points.


ofstream file1 ("X profile y=0.1Ly n=11000.plt");     // output file for the potential profile at timestep n=10.

ofstream file2 ("X profile y=0.1Ly n=12000.plt");     // output file for the potential profile at timestep n=50.

ofstream file3 ("X profile y=0.1Ly n=13000.plt");    // output file for the potential profile at timestep n=50.

ofstream file4 ("X profile y=0.1Ly n=15000.plt");    // output file for the potential profile at timestep n=50.


delx=Lx/Nx; cout<<"delx= "<<delx<<endl;

dely=Ly/Ny; cout<<"dely= "<<dely<<endl; getch();


// ----------------------- Assignment of boundary values ----------------------------------------

for(j=0;j<=Ny;j++){ Phi[0][j]=a4; Phi[Nx][j]=a2;   }

for(i=0;i<=Nx;i++){ Phi[i][0]=a1; Phi[i][Ny]=a3;  }


for(i=0;i<=Nx;i++){ Phinew[i][0]=a1; Phinew[i][Ny]=a3; }

for(j=0;j<=Ny;j++){ Phinew[0][j]=a4; Phinew[Nx][j]=a2; }


// ----------------------- initial potential assignement ----------------------------------------

for(j=1;j<Ny;j++){ for(i=1;i<Nx;i++){ Phi[i][j]=0.1*a3; }}


for (n=1;n<=T;n++){ if (fmod(n,100)==0) cout<<"n= "<<n<<endl; //getch();

for(j=1;j<Ny;j++){
for(i=1;i<Nx;i++){

Phinew[i][j]=0.25*(Phi[i+1][j] + Phinew[i-1][j] + Phi[i][j+1] + Phinew[i][j-1]);

//cout<<"Phinew["<<i<<"]["<<j<<"]= "<<Phinew[i][j]<<endl;

}}

if (n==11000)   { for(i=1;i<=Nx;i++){ file1<<i*delx<<" "<<Phinew[i][Ny/10]<<endl; }}

if (n==12000)   { for(i=1;i<=Nx;i++){ file2<<i*delx<<" "<<Phinew[i][Ny/10]<<endl; }}

if (n==13000)   { for(i=1;i<=Nx;i++){ file3<<i*delx<<" "<<Phinew[i][Ny/10]<<endl; }}

if (n==15000)   { for(i=1;i<=Nx;i++){ file4<<i*delx<<" "<<Phinew[i][Ny/10]<<endl; }}

for(j=1;j<Ny;j++){
for(i=1;i<Nx;i++){

Phi[i][j]=Phinew[i][j];

}}


}

cout<<"end." <<endl; getch();

}